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Abstract 

In this paper, we discuss the incorporation of dynamic subgrid scale (SGS) models in 
the lattice-Boltzmann method (LBM) for large-eddy simulation (LES) of turbulent 
flows. The use of a dynamic procedure, which involves sampling or test-filtering of 
super-grid turbulence dynamics and subsequent use of scale-invariance for two lev- 
els, circumvents the need for empiricism in determining the magnitude of the model 
coefficient of the SGS models. We employ the multiple relaxation times (MRT) 
formulation of LBM with a forcing term, which has improved physical fidelity and 
numerical stability achieved by proper separation of relaxation time scales of hydro- 
dynamic and non-hydrodynamic modes, for simulation of the grid-filtered dynam- 
ics of large-eddies. The dynamic procedure is illustrated for use with the common 
Smagorinsky eddy-viscosity SGS model, and incorporated in the LBM kinetic ap- 
proach through effective relaxation time scales. The strain rate tensor in the SGS 
model is locally computed by means of non-equilibrium moments of the MRT-LBM. 
We also discuss proper sampling techniques or test-filters that facilitate implemen- 
tation of dynamic models in the LBM. For accommodating variable resolutions, we 
employ conservative, locally refined grids in this framework. As examples, we con- 
sider the canonical anisotropic and inhomogeneous turbulent flow problem, i.e. fully 
developed turbulent channel flow at two different shear Reynolds numbers Re* of 
180 and 395. The approach is able to automatically and self-consistently compute 
the values of the Smagorinsky coefficient, Cs- In particular, the computed value 
in the outer or bulk flow region, where turbulence is generally more isotropic, is 
about 0.155 (or the model coefficient C = C\ = 0.024) which is in good agreement 
with prior data. It is also shown that the model coefficient becomes smaller and 
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approaches towards zero near walls, reflecting the dampening of turbulent length 
scales near walls. The computed turbulence statistics at these Reynolds numbers 
are also in good agreement with prior data. The paper also discusses a procedure 
for incorporation of more general scale-similarity based SGS stress models. 

Key words: Lattice-Boltzmann Method, Multiple-Relaxation-Time Model, 
Turbulent Flows, Large Eddy Simulation, Dynamic Subgrid Scale Modeling 
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1 Introduction 



Turbulence in fluids is ubiquitous in nature and technological systems and rep- 
resents one of the most challenging aspects in fluid mechanics. The difficulty 
stems from the inherent presence of many scales that are generally insepara- 
ble, among many other factors. Nevertheless, considerable progress has been 
made over the years towards more fundamental physical understanding of 
turbulence phenomena through measurements, statistical phenomenological 
theories, modeling and computation [T|2|3] . 

In the case of their modeling and computation, different approaches, depend- 
ing on the desired degree of detail, have been devised. One approach is to em- 
ploy Reynolds-averaged equations of fluid motion representing the evolution of 
mean turbulence field. In this framework, the averaged turbulent fluctuations 
- the Reynolds stresses - need to be modeled. Various turbulence models 
have been developed for this purpose, including the well-known k — e and 
Reynolds stress transport equations. A major limitation of this approach is 
that Reynolds-averaged models need to represent physics over a wide range of 
scales. While small scale turbulence tends be somewhat more universal, large 
scale turbulent motions are strongly problem dependent. Thus, it is unreal- 
istic to expect the Reynolds-averaged models to accommodate and represent 
large-scale behavior of different classes of turbulent flows in the same manner 
without resorting to considerable empiricism. 

On the other hand, direct numerical simulation (DNS) approach resolves all 
spatio-temporal scales up to dissipation scales without the use of models, and 
can thus predict all possible motions, structural and statistical features due 
to turbulence with high fidelity. Due to its high computational cost, its appli- 
cation is restricted to relatively low Reynolds numbers. On the other hand, in 
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view of this, it is more pragmatic to only compute the eddying motions of flu- 
ids with length scales greater than the grid size and model unresolved subgrid 
scales (SGS), which comprises the large eddy simulation (LES) approach. The 
SGS scales, which mainly account for dissipation of turbulence, are generally 
more isotropic in nature and relatively independent of the resolved large-scale 
fluid dynamics. Thus, LES represents a compromise with reduced empiricism 
in contrast to Reynolds-averaged models, but with reduced computational cost 
in comparison with DNS. 

In LES approach, the Smagorinsky eddy- viscosity model (1963) [I] is one 
of the simplest SGS models. Based on the assumption that the small scales 
are in equilibrium, they dissipate all the energy transferred from large scales 
instantaneously and completely. This dissipation mechanism is modeled by 
an eddy viscosity, which is represented by the product of the square of a 
length scale, which is often the grid size, and the inverse of a time scale, which 
is given by the local fluid shear rate, with a proportionality coefficient or 
"constant" . The values of this proportionality "constant" were first determined 
by Lilly (1966) [5] and Deardorff (1970) [6] by analysis for isotropic turbulence 
and turbulent channel flow, respectively, and one of the first successful LES 
applications is by Schumann (1975) [7]. A more theoretical basis for LES was 
first provided by Leonard (1974) [8] through the introduction of formal filtering 
procedures on fluctuating turbulence fields. 

In wall-bounded flows, turbulence becomes anisotropic close to walls, with its 
length scales becoming progressively smaller and vanishing at the wall. Thus, 
with the use of a finite grid size as the length scale in the SGS model, the 
proportionality coefficient should decrease and approach to zero towards the 
wall to account for near-wall turbulence anisotropy and the model coefficient 
is actually not a "constant". Often, an ad hoc damping function [9] is used 
for this purpose [TO]. Moreover, in transitional flows, it is found the use of 
the "constant" Smagorinsky SGS model in LES can cause excessive damp- 
ing of the resolved structures [11] and can be alleviated by an additional ad 
hoc intermittency function. The need to a priori specify the values of the 
proportionality coefficient and their variation, for example, in the presence of 
turbulence anisotropy or laminar-to-turbulence transition comprises the em- 
pirical elements of the LES approach. 

Significant progress in the LES approach was made by the introduction of 
dynamic procedures for local computation of the coefficient (s) in the SGS 
models by Germano and co-workers [T2fT5] . By obtaining information of the 
turbulence field at scales larger than the grid size - the "test" filter scale, 
and relating them with the resolved turbulence at grid scale through scale- 
invariance, the value of the model coefficient can be computed dynamically 
as the computation progresses. Thus, dynamic procedures avoid the need to 
a priori specify the model coefficient and determine them based on the local 
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fluctuating turbulence field. Moreover, unlike the constant Smagorinsky SGS 
model, it does not rule out the possibility of backscattering of energy from 
smaller scales. It is found that the use of dynamic procedure makes LES 
susceptible to numerical instabilities, which can be substantially alleviated 
by employing a least squares technique [H] in conjunction with the use of 
averaging to compute the coefficient. 

When the dynamic procedure is applied to the Smagorinsky SGS model, it 
is often simply referred to as the dynamic Smagorinsky model (DSM). Other 
SGS models include the use of a scale-similarity stress tensor combined with 
the Smagorinksy model - the mixed SGS model [15j - to provide improved 
predictions for certain classes of inhomogeneous flows. Dynamic procedures 
have been applied with much success to mixed SGS models, resulting in dif- 
ferent variants such as the dynamic mixed model (DMM) [T6] and the dynamic 
two-parameter model (DTM) [TTJ . The remarkable progress made in the devel- 
opment of LES as a technique for turbulence simulation and their applications 
have been discussed in detail in the reviews by Piomelli (1999) [18], Meneveau 
and Katz (2000) [19] and Pope (2004) [20], and in the monographs of Galperin 
and Orszag (1993) [2TJ and Sagaut (2003) [22]. 

In this paper, we propose to incorporate dynamic procedure in the lattice 
Boltzmann method (LBM) for LES. In recent years, the LBM, based on mini- 
mal discrete kinetic models, has emerged as an alternative and accurate com- 
putational approach for fluid mechanics problems [2"3"f2"4] . It involves the solu- 
tion of the lattice-Boltzmann equation (LBE) that represents the changes in 
the evolution of the distribution of particle populations due to their advection, 
represented as a free-flight process, and collision, described as a relaxation pro- 
cess, on a lattice. When the lattice, which represents the discrete directions 
for particle propagation, satisfies sufficient rotational symmetries, the averaged 
LBE asymptotically recovers the weakly compressible Navier-Stokes equations 
(NSE). 

Though it evolved as a computationally efficient form of lattice gas cellular 
automata [25] , it was well established about a decade ago that the LBE is ac- 
tually a much simplified form of the Boltzmann equation [26,27] . As a result, 
several previous results in discrete kinetic theory could be directly applied to 
the LBE. This lead to, for example, improved physical modeling in various 
situations, such as multiphase flows [25129] and multicomponent flows |30j . 
and in an asymptotic theory suitable for rigorous numerical analysis [3T] • As a 
result of features of the stream-and-collide procedure of the LBE such as the 
algorithmic simplicity, amenability to parallelization with near-linear scalabil- 
ity, and its ability to represent complex boundary conditions and incorporate 
physical models more naturally, it has rapidly found a wide range of applica- 
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Discrete kinetic theory, and in particular the LBM, has recently been em- 
ployed to a priori derive turbulence models based on a mean-field approxi- 
mation [36,3 71138] and lattice kinetic turbulence models have found much suc- 
cess |39j. It has now been well established that LBM is a reliable and accurate 
method for direct numerical simulation (DNS) of various benchmark turbu- 
lent flow problems [4T)ll41ll42f43|l44j . On the other hand, in earlier work, the use 
of LBM as a LES tool with "constant" Smagorinsky model to represent SGS 
effects was also proposed [^511^6] . which has more recently found applications 
for flows in different configurations and physical conditions |4"7I3 8,49,50]. 

A commonly used form of the LBM employs a single relaxation time (SRT) 
model [51] to represent the effect of particle collisions, in which particle dis- 
tributions relax to their local equilibrium at a rate determined by a single 
parameter [52"f53"] . On the other hand, an equivalent representation of distri- 
bution functions is in terms of their moments, such as various hydrodynamic 
fields including density, mass flux, and stress tensor. The relaxation process 
due to collisions can more naturally be described in terms of a space spanned 
by such moments, which can, in general relax at different rates. This forms the 
basis of the generalized lattice-Boltzmann equation (GLBE) based on multiple 
relaxation times (MRT) [54"j55ll56j . 

By carefully separating the time scales of various hydrodynamic and kinetic 
modes through a linear stability analysis, the numerical stability of the GLBE 
or MRT-LBE can be significantly improved when compared with the SRT- 
LBE, particularly for more demanding problems at high Reynolds numbers [55] . 
MRT-LBE has also been extended for multiphase flows with superior stabil- 
ity characteristics [57j58ll59ll6"0] , and, more recently, for magnetohydrodynamic 
problems [SI] . It has also been used for LES of a class of turbulent flows [4"8||50] . 
Recently, we have extended the GLBE with forcing terms for wall-bounded 
turbulent flows [62J, which is a generalization of the approach introduced ear- 
lier by Yu et al. [50J. These forcing terms, determined to avoid discrete lattice 
artifacts, represent the effect of external forces and are considered in nat- 
ural moment space of GLBE. In this work, we represented the SGS effects 
through the relaxation times of the hydrodynamic moments using the "con- 
stant" Smagorinsky eddy- viscosity model [I], which is modified by the van 
Driest damping function [9] in the near wall region. 

One of the main objectives of this paper is to develop and implement a frame- 
work for application of dynamic procedure in the lattice Boltzmann method 
(LBM) for LES of inhomogeneous and anisotropic turbulent flows. In view of 
the previous discussions, we plan to employ GLBE based on multiple relax- 
ation times as the basis for this purpose. We will present the incorporation of 
the dynamic procedure for the Smagorinsky SGS model, i.e. DSM, in LBM. 
Particular attention is paid to performing discrete filtering operations at su- 
per lattice grid or test-filter levels, which will be discussed. To accommodate 
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efficient simulation of varying turbulent length scales in the near-wall region, 
variable resolution strategies are required. In this regard, we employ conser- 
vative, locally refined multiblock grids that allow variable resolution. The ap- 
proach will be tested for the fully-developed turbulent channel flow problem 
at two different shear Reynolds numbers i?e* of 180 and 395, by comparing 
the turbulent statistics with prior data based on direct numerical simulations 
(DNS) and measurements. The higher Reynolds number case, to the best of 
our knowledge, is considered for the first time in the LBM for LES of tur- 
bulent channel flows. The incorporation of DSM in the GLBE is a first step 
in the application of dynamic procedure. More general SGS models, involving 
scale-similarity stresses, such as DMM and DTM can also be introduced in 
the GLBE, which will also be briefly described in an appendix of this paper, 
and their implementation will be a subject of future work. 

The organization of this paper is as follows: In Sec. [2J we will discuss the 
generalized lattice-Boltzmann equation (GLBE) with forcing term. The use of 
conservative locally refined multiblock grids will be presented in Sec. [3j Sec- 
tion H] will discuss the dynamic SGS modeling with DSM in GLBE framework. 
Then, Sec. [5] will present the LES results for the canonical turbulent channel 
flow problem at two different Reynolds numbers. Finally, the summary and 
conclusions of this paper will be presented in Sec. El 



2 Generalized Lattice Boltzmann Equation with Forcing Term 



We shall now discuss the computational procedure based on generalized lattice- 
Boltzmann equation with a forcing term, which is supplemented by Smagorin- 
sky SGS model with a variable coefficient, obtained by means of a dynamic 
procedure (see next section). For brevity, we will present only the major ele- 
ments of the approach, while the details can be found in Refs. [56.50.58.62j. 

The lattice-Boltzmann method computes the evolution of distribution func- 
tions as they move and collide on a lattice grid. The collision process consider 
their relaxation to their local equilibrium values, and the streaming process 
describes their movement along the characteristics directions given by a dis- 
crete particle velocity space represented by a lattice. Figure [1] represents the 
three-dimensional, nineteen particle velocity (D3Q19) lattice model employed 
in this paper. The particle velocity corresponding to this model may be 
written as: 



(0,0,0) a = 

(±1,0,0),(0,±1,0),(0,0,±1) « = l,---,6 
(±1, ±1, 0), (±1, 0, ±1), (0, ±1, ±1) a = 7, ■ ■ ■ , 18 
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Fig. 1. Schematic illustration of the three-dimensional, nineteen velocity (D3Q19) 
model. 



The GLBE computes collision in moment space, while the streaming process 
is performed in the usual particle velocity space [56]. The GLBE with forcing 
term [52] also computes the forcing term, which represents the effect of exter- 
nal forces as a second-order accurate time-discretization, in moment space. We 
use the following notation in our description of the procedure below: In particle 
velocity space, the local distribution function f , its local equilibrium distribu- 
tion f eq , and the source terms due to external forces S may be written as the 
following column vectors: f = [/„, fa fa . . . , fa}\ f e? = [/o 9 , f{ q , f% q , • • • , ffgf, 
and S = [So, Si, S 2 , ... , Si 8 ] ■ Here, the superscript f represents the transpose 
operator. 

The moments f are related to the distribution function f through the rela- 
tion f = Tf where T is the transformation matrix. Here, and in the follow- 
ing, the "hat" represents the moment space. The transformation matrix T is 
constructed such that the collision matrix in moment space A is a diagonal 
matrix through A = TAT -1 , where A is the collision matrix in particle ve- 
locity space. The elements of T are obtained in a suitable orthogonal basis 
as combinations of monomials of the Cartesian components of the particle 
velocity through the standard Gram-Schmidt procedure, which are pro- 
vided by d'Humieres et al. (2002) [56]. Similarly, the equilibrium moments 
and the source terms in moment space may be obtained through the trans- 
formation i eq = Ti eq , S = TS. The components of moment-projections of 

these quantities are: f = f , fa fa . . . , fa , i eq = f^, f* q , f 2 q , . . . , f eq 



is 



So, Si, S 2 , ■ ■ ■ , Sis ■ These are provided in Appendix \M 



t 



and S 

The solution of the GLBE with forcing term can be written in terms of the 



7 



following "effective" collision and streaming steps, respectively: 



f(l?,t) =f(z>,t) +mC#,t), (2) 



and 

f a (l? + ~e> a 8 t ,t + 8t) = f a (!?,t), (3) 

where the distribution function f = {/ Q } Q =o,i,...,i8 is updated due to "effective" 
collision resulting in the post-collision distribution function f = {/ a }a=o,i,...,i8 
before being shifted along the characteristic directions during streaming step. 
vj represents the change in distribution function due to collisions as a relax- 
ation process and external forces, and following Premnath et al. (2008) [52] it 
can written as 



-A (f-f '')_> +(x--fi)s,-=>~ 



(4) 



where X is the identity matrix and A = diag(s , si, . . . , Sig) is the diagonal 
collision matrix in moment space. Now, some of the relaxation times s a in 
the collision matrix, i.e. those corresponding to hydrodynamic modes can be 
related to the transport coefficients and modulated by eddy-viscosity due to 
SGS model (discussed below) as follows [56|50f62j : sj -1 = |C + |> where ( 
is the molecular bulk viscosity, and sg = su = S13 = su = S15 = s„, where 
s^ 1 = 3z/+| = 3(z/o + ^) + |j with v being molecular shear viscosity and v t the 
eddy-viscosity determined from the dynamic SGS model discussed in the next 
section. The rest of the relaxation parameters, i.e. for the kinetic modes, can be 
chosen through a von Neumann stability analysis of the linearized GLBE |56j : 
s\ = 1.19, S2 = S10 = S12 = 1-4, S4 = sq = s$ = 1.2 and s\q = Syj = si& = 1.98. 

Once the distribution function is known, the hydrodynamic fields, i.e., the 
density p, velocity it and pressure p can be obtained as follows: 

18 18 1 _ 

P=^2fa, T = P = X! foT&a + ~F8t, V = C 2 s p } (5) 
a=0 a=0 1 



where, c s = cj v3 with c = 5 X / 5 t being the particle speed, and 8 X and S t are the 
lattice spacing and time step, respectively. It can be readily shown that when a 
multiscale analysis based on the Chapman- Enskog expansion [63.58J is applied 
to the GLBE with relaxation times scales augmented by an eddy-viscosity, it 
recovers the "grid filtered" weakly compressible Navier- Stokes equations with 
external forces. That is, the density, momentum and pressure obtained from 
Eq. (J5J) are grid-filtered quantities: p — > p,plt — > pu and p — > p, where 
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the 'overbar' represents grid-filter. It may be noted that the macrodynamical 
equations can alternatively be derived through an asymptotic analysis under 
a diffusive scaling [51] . 

The computational procedure for the solution of the GLBE with forcing term 
is optimized by fully exploiting the special properties of the transformation 
matrix T: these include its orthogonality, entries with many zero elements, 
and entries with many common elements that are integers, which are used 
to form the most common sub-expressions for transformation between spaces 
in avoiding direct matrix multiplications [55]. For details, we refer the reader 
to Ref. [5T)1IB"2] . As a result of such optimization, the additional computational 
overhead when GLBE is used in lieu of the popular SRT-LBE is small, typically 
between 15% — 30%, but with much enhanced numerical stability that allows 
maintaining solution fidelity on coarser grids and also in simulating flows at 
higher Reynolds numbers. 



3 Conservative Local Grid Refinement Based on Multiblock Ap- 
proach 

Closer to a wall, length scales are very small, requiring a fine grid to adequately 
resolve turbulent structures. Use of grid fine enough to resolve the wall regions 
throughout the domain can entail significant computational cost, and this can 
be mitigated by introducing coarser grids farther from the wall, where turbu- 
lent length scales are larger. One approach is to consider using continuously 
varying grid resolutions, using a interpolated-supplemented LBM [61] that ef- 
fectively decouples particle velocity space represented by the lattice and the 
computational grid. However, it is well known that interpolation could intro- 
duce significant numerical dissipation, see for e.g. [55], which could severely 
affect the accuracy of solutions involving turbulent fluctuations. Thus, we 
consider locally embedded grid refinement approaches, and in particular their 
conservative versions [65.66J that preserve mass and momentum conservation. 
Similar zonal embedded approaches have been successfully employed in com- 
putational approaches based on the solution of filtered NSE for LES of turbu- 
lent flows [67] . 

Figure [2] shows a schematic of such a multiblock approach in which a fine 
cubic lattice grid is used close to the bottom wall surfaces and a coarser one, 
again cubic in shape, farther out. Similarly, finer grids are employed near top 
free-slip surfaces. In order to facilitate the exchange of information at the 
interface between the grids, the spacing of the nodes changes by an integer 
factor, in this case two. As well as using different grid sizes, the two regions 
use different time steps (time step being proportional to grid size), and the 
computational cost required per unit volume is thus reduced by a factor of 
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Fig. 2. Schematic illustration of conservative, locally refined multiblock gridding 
strategy for LES of turbulent channel flow using GLBE with dynamic SGS model. 

16 in the coarse grid. Figure [2] shows a staggered grid arrangement, in which 
nodes on the fine and coarse sides of the interface are arranged in a manner 
that facilitates the imposition of mass and momentum conservation. Different 
blocks communicate with each other through the Coalesce and the Explode 
steps, in addition to the standard stream- and- collide procedure. The details 
are provided in Refs. [65.66J and we very briefly present the essential elements 
in the following. 

The Coalesce procedure involves summing the particle populations on the 
fine nodes to provide new incoming particle populations for the corresponding 
coarse nodes. Similarly, the Explode step involves redistributing the popula- 
tions on the coarse node to the surrounding fine nodes. These grid - com- 
municating steps used in the multiblock approach presented in Ref. [65J were 
incorporated in GLBE framework in this work. It may be noted that other 
approaches to local grid refinement exist in the literature. In particular, asymp- 
totic analysis of the LBE shows that it converges to the incompressible NSE 
when the ratio of the square of time step and square of the grid spacing is 
a constant [21]. Local grid refinement based on this consideration need to 
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adapt the time step to vary as the square of the grid spacing. However, in this 
work we use the grid refinement variant that considers linear variation of the 
time step with the grid spacing for reasons of relatively lower computational 
cost and the fact that it is already more developed and tested for different 
problems [65|66] . 




4 Subgrid Scale Eddy- Viscosity Modeling using Dynamic Proce- 



The main idea behind the application of dynamic procedure to SGS mod- 
els [12] is that the local value of the model coefficient of the eddy viscosity, 
representing eddy-viscosity type subgrid scale effects, can be obtained from 
sampling the smallest super-grid or resolved scales, which are generally re- 
ferred to as the test-filtered scales, and assuming scale-invariance at these 
two levels. If A is the width of the grid filter, which in the GLBE becomes 
A = 8 X , where 5 X is the lattice grid spacing, the flow information is sampled at 
a larger scale A, the test-filter width, i.e. A > A, and generally A/ A = 2 |12j . 
The notation adapted here, as before, and in the following is that 'bar' refers 
to grid-filtered values and a 'tilde' refers to test-filtered values. The effect of 
subgrid scales is parameterized by an eddy viscosity relation [I] 



where \S\ is the strain rate, given by |>S| = y2SijSij and C is the model 

coefficient. The grid-filtered strain rate \S\ can be obtained by evaluating the 
corresponding strain rate tensor from the non-equilibrium moments as shown 
in Appendix [A] The model coefficient C is obtained as follows. 

The anisotropic part of the SGS stress at grid- filter scale Ty (r^ = UiUj —HiU]) 
and that at the test-filter scale (T^- = UiUj — u^u]), where UjU] and UjU] 
are unknowns, are modeled in terms of eddy viscosity and strain rates at 
corresponding scales, respectively, or by invoking scale-invariance, as 




dure 



u t = CA 2 \S\ 



(6) 




Ti 




(7) 



and 



where Sy and S refers to the strain rate tensors at the grid- and test- filter 

levels, respectively. The respective magnitude of strain rates \S\ and \S\ may 
be written as: 

\S\ = ^SijSij, Sij = 1/2 (diUj + djUi), 
\S\ = \j2SijSij, = 1/2 (diUj + djU^j. 

It may be noted that the grid-filtered strain rates can be obtained locally 
from the non-equilibrium moments, as mentioned above; the test-filtered strain 
rates can be determined directly from the gradient of test-filtered velocity 
field, which can, in turn, be obtained by explicit application of test-filter on 
the grid-filtered velocity field. The expressions for the test-filter that can be 
conveniently used in the GLBE framework will be discussed later. 

The unknown SGS stress at each filter level can be related by the Germano 
identity pj] 

Lij UiUj UiUj ^ij Tij , (9) 



where Ly is the resolved turbulent stress, and upon substituting Eqs. (j7j) 
and (|SJ) in Eq. Q), we get an expression between Ly and other known variables 
with the model coefficient C being the only unknown 

Lij - -fL kk = -2CM tv (10) 



where 



Mij = A \S\S tJ - A 2 \S\S lr (11) 



We note that the computation of the first term in the My tensor involves 
explicit test-filtering and finite-differencing, while its second term can be ob- 
tained locally from non-equilibrium moments, as discussed above. Since, the 
tensor expressions, i.e. Eq. (fTOj) and (iTTj) leads to five independent equations 
containing C, it can be obtained by a least square minimization approach 
as [Hj: 

c 1 (LkiM kl ) 

2(M i3 MijY 1 1 



where the usual summation convention of the repeated indices is assumed and 
< ■ > represents spatial averaging in homogeneous directions or time averaging 
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or both, depending on the problem, so as to stabilize the computations. Also 
clipping is to be done when u t < 0, i.e. set v t = in such cases. Once the 
eddy viscosity in Eq. (jBJ) has been determined from Eqs. (19|)- (IT2|) it is added 
to the molecular viscosity u , obtained from the statement of the problem, 
to yield the total viscosity u, i.e. v = v$ + u t , which can then be used to 
compute the "effective" relaxation times for hydrodynamic moments for the 
GLBE discussed in Sec. |2j 



4-1 Test- Filters 



Several approaches exist to obtain the test-filtered quantities from the grid- 
filtered variables that can, in turn, be directly obtained from the solution of 
the GLBE. The approach considered in this paper is the discrete trapezoidal 
filter [TB,68j. Discrete trapezoidal filters are employed in this work as they 
naturally accommodate with the underlying cubic lattice grid structure and 
thereby allowing an efficient implementation. They involve applying the one- 
dimensional filter successively in each of the three coordinate directions as 
follows: 

&*i,j,k) = V 4 {<f>(i+i,j,k) + 2 0(i,j,k) + 0(i-ij,fc)) > 

= V 4 (%,k+i) + 2 ^*tj,k) + 0*tj,k-i)) , 

to obtain the test-filtered value (f> from at (i,j,k). Near walls, along the 
coordinate direction k, the test-filtered values at the top and bottom can be ob- 
tained as j = 1/2 (0(* i)fc) + 0(* i>fe+ i)) and0 (ij - fc) = 1/2 » 
respectively. 

Other types of test-filters specialized for lattice stencils that can also maintain 
isotropy may also be possible, but presumably at the expense of additional 
computational overhead. The implementation of test-filters in the multiblock 
approach (see Sec. [3]) leads to additional challenges in how best to deal with 
performing such filtering operations at the grid boundaries. In order to cal- 
culate the model coefficient C, it is necessary to use information from two 
grid spacings away; thus to compute C at the first node on the coarse grid, 
two ghost nodes are introduced. Velocities at these ghost nodes are obtained 
by averaging the values at the surrounding fine nodes. Note that since the 
collision step is not performed on the two fine layers adjacent to the grid 
boundary, it is not necessary to introduce additional ghost nodes for the fine 
grid. For an efficient implementation, optimization strategies are considered 
in the incorporation of the dynamic procedure in the multiblock GLBE. 
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It may be noted that the dynamic procedure with other SGS models, such 
as the dynamic mixed model (DMM) [16] and the dynamic two-parameter 
model (DTM) [T7] can be incorporated in the GLBE in a similar manner. 
The details of the incorporation of these models in the GLBE are presented 
in Appendix [B] 



5 Results and Discussion 

We will now present investigations of the LBM using the GLBE formulation 
with dynamic SGS modeling for LES of wall-bounded turbulent flows. In ear- 
lier studies, the LBM has been validated as a direct numerical simulations 
(DNS) approach for such problems [4"T1l4"l] . In the present study, we will assess 
the performance of GLBE as a LES approach that incorporates SGS effects 
using a dynamic procedure to compute the model coefficient, thereby avoid- 
ing empiricism, on a relatively coarse grid with multiblock approach, while 
maintaining necessary near-wall resolutions. 

5.1 Turbulent Channel Flow at Re* = 180 

First, we consider fully-developed turbulent flow in an open channel with a 
shear Reynolds number of Re* = u*H/ 1/ = 180, where u* is the shear velocity, 
vq is the molecular kinematic viscosity and H is the channel height. The shear 
velocity is related to the wall stress r w by u* = Jr w /p. It represents an 
important canonical problem to test turbulence models and computational 
techniques, for which detailed direct numerical simulations (DNS) data are 
available for comparison [691170] . 

5.1.1 Computational Conditions 

No-slip boundary is considered at the bottom, implemented using half-way 
link bounce back approach [71] . free-slip at the top, by means of specular 
reflection of distribution functions |24J, and periodic boundary conditions are 
applied in the streamwise, x and spanwise, y directions. The computational 
domain is chosen with appropriate aspect ratios, viz., 5.7 H and 2.85H in the 
streamwise and spanwise directions, respectively. With these aspect ratios, a 
sufficient number of wall-layer streaks are accommodated and end effects of 
two-point correlations are excluded, i.e. the two-point velocity correlations in 
solutions are required to decay nearly to zero within half the domain [72]. 

For wall-bounded turbulent flows, it is important to adequately resolve the 
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near-wall, small-scale turbulent structures. Grid spacing for such problems are 
generally referred in terms of wall units (designated with a "+" superscript) as 
A + = A/8 V , where A is the wall spacing and 5 U = vo/u* is the characteristic 
viscous length scale. When the computations resolve the local dissipative or 
Kolmogorov length scale rj = (z/p/e) 1 / 4 , i.e. A+ w < 0(r/ + ) [72], it can satisfy 
the near-wall resolution requirements. In particular, it is generally recognized 
that 1.5?7 + — 2. Or/ 4 " represents the upper limit of grid-spacing, above which the 
small scale turbulent motions in bounded flows are not well resolved. It can 
be shown by simple arguments that rj + « 1.5 — 2.0 at the wall and that i] + 
increases with increasing distance from the wall [2J. 

To account for resolution requirements, we discretize the computational do- 
main in terms of two blocks with different resolutions, viz., a fine block near 
the wall and a coarse block, whose grid spacings are twice that of the fine- 
block in each direction, in the bulk bounded by free-slip surface at the top. 
For the fine-block, we used A^ ine = 3.95 in wall units. Due to the use of 
link-bounce back method for implementation of wall boundary condition, the 
first lattice node is located at a distance of A+ w = A + /2, which in our case 
is 1.98. Hence, the computations are expected to fairly resolve the small-scale 
turbulent structures. For the coarse-block, we used A+ arse = 2 x A^ ine = 7.90. 
The computational domain is discretized by 256 x 128 x 12 grid nodes and 
128 x 64 x 19 grid nodes in the fine and coarse blocks, respectively. 

The initial mean velocity is specified to satisfy the I /7 th power law [2], while 
initial perturbations satisfying divergence free velocity field [70]. The density 
field is taken to be p = po — 1-0 f° r the entire domain. The precise form of the 
initial fields may not affect the turbulence statistics, but can have significant 
influence on the convergence of the solution to statistically steady state. In 
particular, the above choice of initial fields would enable a rapid convergence 
to the statistically steady state solution of the fluctuating fields obtained by 
GLBE with forcing term. We employed the consistent initialization proce- 
dure for the distribution functions or moments [73], that properly initializes 
non-equilibrium moments of the GLBE in the presence of non-uniform hydro- 
dynamic fields, such as those mentioned above. 



5.1.2 Computational Procedure 

Using ~P = — ^x = ^x = ^jfx as the driving force, LES using the GLBE 
is performed. The collision step, including the forcing term, is computed in 
moment space, while the streaming step is carried out in its natural velocity 
space, and the result of these two steps provides the grid-filtered hydrody- 
namic fields. The grid-filtered strain rate tensor Stj is obtained locally from 
non-equilibrium moments as shown in Appendix [A], which was also found to 
provide improved numerical stability on coarser grids when used in lieu of 
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finite-differencing of the velocity fields. The relaxation process in the collision 
step is carried out by computing the "effective" relaxation times for hydrody- 
namic moments obtained from a variable coefficient Smagorinsky SGS eddy- 
viscosity model, whose values are determined by a dynamic procedure. The 
steps to accomplish this is briefly summarized as follows. 

Based on these grid- filtered quantities and considering A/A = 2 where A = 5 X , 
discrete trapezoidal test-filters, as discussed in Sec. I4.1[ are applied to obtain 

the following quantities: Ui, UiUj, Sij, \S\ and l^lSy. Once they are known, 
tensors such as and used in the dynamic procedure (Eqs. ([9]and[TT]), re- 
spectively, see Sec. HI) are computed. The model coefficient C is, then, obtained 
by averaging the contracted tensor products, as shown in Eq. ( |T2|) . along the 
homogeneous x and y directions, i.e. horizontal planes, to maintain numerical 
stability. Moreover, near free-slip surfaces, we specify that there is no change 
in the model coefficient around the top surface to ensure stability. Other SGS 
models, such as the DTM whose implementation details in GLBE framework 
are provided in the Appendix [HI could be used to relax this assumption, which 
is a subject for a future work. Finally, the SGS eddy-viscosity is obtained from 
Eq. ([6]) to provide the "effective" relaxation times for the collision step. The 
GLBE computations are carried out until stationary turbulence statistics are 
obtained, as measured by the invariant Reynolds stresses profiles. This initial 
run was carried out for a duration of 50T*, where T* = H/u* is the charac- 
teristic time scale. The averaging of various flow quantities was carried out in 
time as well as in space in the homogeneous directions by an additional run 
for a period of 12T*. 



5.1.3 Results 

Figure [3] shows the variation of the SGS model coefficient Cd yn = C along 
the wall normal direction, where the distance is scaled in the wall coordinates 
Z + , where Z + = Z/5 U and where 5 U is the viscous length scale defined earlier. 
The model coefficient is seen to be nearly a constant in the bulk or outer flow 
region, where turbulence tends towards becoming statistically more homoge- 
neous and isotropic. In particular, the value of Cd yn in these regions is found 
to be about 0.024. This is excellent agreement with Germano et al. (1991) [12] 
based on the solution of filtered Navier-Stokes equations, who found a value 
of 0.023. The eddy- viscosity given in Eq. ([6]), is also written in terms of the 
Smagorinsky "constant" C$ as v t = Cf A \S\. Thus, C,s = VC and we find 
Cs = 0.155 in the bulk region, which is within the range of 0.10 — 0.16 reported 
for various classes of flows [21] . Near Z + = 40, the grids transition from coarse 
to fine block, where some oscillations in the values of C are noticed, which 
could be removed by an additional averaging. Closer to the wall (Z + < 30), 
in the inner-layer region, the model coefficient is found to decrease monoton- 
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Fig. 3. Computed variation of model coefficient for fully-developed turbulent channel 
flow with a top slip surface at a shear Reynolds number of i?e* = 180 using a 
dynamic procedure with multiblock GLBE. 



ically with decreasing distance from the wall Z + and approaching zero very 
close to the wall. This trend is consistent with near-wall turbulence physics, 
according to which turbulence is statistically anisotropic and its scales become 
smaller towards the wall. Thus, the GLBE using dynamic Smagorinsky SGS 
model is not only able to self-consistently predict the bulk value of the model 
coefficient automatically and accurately, but also its variation near the wall 
without resorting to any empirical approach such as the van Driest damp- 
ing function [9]. Such predictive capabilities are among the major assets of 
dynamic SGS modeling using GLBE. 

Figure H] shows the computed mean velocity profile, normalized by the shear 
velocity function of the distance from the wall given in wall units, i.e. 

Z + . Also plotted are the DNS data by Kim, Moin and Moser (1987) based 
on a spectral method [69]. The computed velocity profile follows fairly closely 
with DNS data, with about 5% difference. Such differences are characteristic of 
LES, which employ relatively coarser grids than DNS, and which also generally 
depends on the numerical dissipation of the computational approach for LES 
(see e.g., Ref. [73175] ). 

The Reynolds stress, normalized by the wall-shear stress, is presented in Fig. [3J 
in semi-log scale and compared with the DNS data of Ref. JHl] obtained 
from the direct solution of incompressible Navier-Stokes equations (NSE). 
It is found that, in general, the computed Reynolds stress is in reasonably 
good agreement with DNS. The computed values are somewhat under pre- 
dicted near Z + = 40, where the transition between coarse to fine blocks of 
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Fig. 4. Mean velocity in fully-developed turbulent channel flow with a top slip 
surface at a shear Reynolds number of Re* = 180 using a dynamic procedure with 
multiblock GLBE. Comparison with DNS of Kim, Moin and Moser (1987). 




Fig. 5. Mean Reynolds stress in fully-developed turbulent channel flow with a top 
slip surface at a shear Reynolds number of Re* = 180 using a dynamic procedure 
with multiblock GLBE. Comparison with DNS of Kim, Moin and Moser (1987). 

grid nodes occur. It is well known that where different types of grids interface, 
certain amount of numerical effects are inevitable, as was found in certain 
flow problems [66]. Nevertheless, it is clear that the computations are able to 
predict the Reynolds stress variation both qualitatively and quantitatively in 
a reasonably accurate manner. 
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Fig. 6. Root-mean-square(rms) streamwise velocity fluctuations in fully-developed 
turbulent channel flow with a top slip surface at a shear Reynolds number of 
Re* = 180 using a dynamic procedure with multiblock GLBE. Comparison with 
DNS of Kim, Moin and Moser (1987) and measurements of Kreplin and Eckelmann 
(1979). 

Second-order statistics, such as turbulence intensities are important measures 
of turbulent activity and let us now assess how the LES computations using 
the GLBE with a dynamic SGS model are able to predict such quantities in 
the near-wall region. Figures El [7J and El show comparisons of the components 
of the root-mean-square (rms) streamwise, spanwise and wall-normal veloc- 
ity fluctuations, respectively, computed using the GLBE with data from DNS 
based on the NSE [69] and experimental measurements [76]. Again, the com- 
puted results are in good agreement with prior data. It is interesting to note 
among the three components of turbulence intensities, only the wall normal 
component is the most sensitive to the numerical effects due to coarse-to-fine 
grid interfaces, which are in any case reasonably small. 



5.2 Turbulent Channel Flow at Re* = 395 



It is important to know whether the approach discussed in this paper is able 
to scale well and reproduce turbulence statistics at higher Reynolds num- 
bers. In this regard, we assess the LES approach based on GLBE using dy- 
namic SGS modeling at a higher shear Reynolds number of 395, for which 
DNS data is available from a more recent work by Moser, Kim and Mansour 
(1999) [77] . To the best of authors' knowledge, LES using LBM, even with 
constant Smagorinsky SGS model, at such a higher Reynolds number and 
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Fig. 7. Root-mean-square (rms) spanwise velocity fluctuations in fully-developed tur- 
bulent channel flow with a top slip surface at a shear Reynolds number of i?e* = 180 
using a dynamic procedure with multiblock GLBE. Comparison with DNS of Kim, 
Moin and Moser (1987) and measurements of Kreplin and Eckelmann (1979). 
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Fig. 8. Root-mean-square(rms) wall normal velocity fluctuations in fully-developed 
turbulent channel flow with a top slip surface at a shear Reynolds number of 
-Re* = 180 using a dynamic procedure with multiblock GLBE. Comparison with 
DNS of Kim, Moin and Moser (1987) and measurements of Kreplin and Eckelmann 
(1979). 
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Fig. 9. Computed variation of model coefficient for fully-developed turbulent channel 
flow with a top slip surface at a shear Reynolds number of Re* = 395 using a 
dynamic procedure with multiblock GLBE. 



their comparison with the DNS has not previously been reported in the liter- 
ature. In order to properly sample sufficient number of wall-layer streaks and 
exclude end effects of two-point velocity correlations, we choose a larger com- 
putational domain than in the previous case, with a size of 5.8H x 2.9H x H. 
The computational domain is discretized in terms of four blocks, with the 
finest grid block near walls, with a spacing of A + = 3.95. The first grid node 
is located at A+ TO = A + /2 = 1.96 due to the use of half-way link bounce back, 
and hence can resolve the smaller near-wall turbulence structures properly. 
The four blocks are resolved by means of the following number of grid nodes: 
576 x 288 x 12, 288 x 144 x 12, 144 x 72 x 15 and 288 x 144 x 8, and the GLBE 
computations are carried out till stationary turbulence are observed, followed 
by an additional run for a duration of 13H/u* for obtaining turbulence statis- 
tics. 



Figure [9] shows the computed variation of the model coefficient C^ yn = C as a 
function of the distance from the wall Z + for Re* = 395. Again, the dynamic 
LES approach is able to predict the value of the model coefficient in the bulk 
or outer flow region accurately. Moreover, it is also able to reproduce the 
near- wall trend of C& yn well, viz., its reduction as the distance from the wall 
is decreased. It is also noticed that around regions, where there is grid block 
transition, some oscillations in C& yn are present. The magnitude of these oscil- 
lations appear to progressively increases with grid transitions involving coarser 
grid blocks, with the maximum occurring at around Z + = 120. However, they 
are generally small and do not seem to affect the turbulent statistics as shown 
below. 
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Fig. 10. Mean velocity in fully-developed turbulent channel flow with a top slip 
surface at a shear Reynolds number of Re* = 180 using a dynamic procedure with 
multiblock GLBE. Comparison with DNS of Moser, Kim and Mansour (1999). 



GLBE 
DNS - 



MKM (1999) 



0.8 




0.6 



0.4 




400 



Fig. 11. Mean Reynolds stress in fully-developed turbulent channel flow with a top 
slip surface at a shear Reynolds number of Re* = 395 using a dynamic procedure 
with multiblock GLBE. Comparison with DNS of Moser, Kim and Mansour (1999). 

Figure [TU] shows the computed mean velocity profile, normalized by the shear 
velocity u*, while the Reynolds stress, normalized by the wall shear stress, for 
Re* = 395 is presented in Fig. [TTJ The DNS data of Moser, Kim and Mansour 
(1999) [77J are also plotted in these figures for comparison. Notice that the 
peak Reynolds stress values are higher for the Re* = 395 case as compared 
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Fig. 12. Root-mean-square (rms) streamwise velocity fluctuations in fully-developed 
turbulent channel flow with a top slip surface at a shear Reynolds number of 
Re* = 395 using a dynamic procedure with multiblock GLBE. Comparison with 
DNS of Moser, Kim and Mansour (1999). 

to the Re* = 180, reflecting greater turbulent activity at the higher Reynolds 
number. It is found that the GLBE LES approach using dynamic SGS model 
is able to reproduce these quantities reasonably well. 

Let us now present the root-mean-square (rms) velocity fluctuations in the 
streamwise, spanwise and wall normal directions in Figs. [12} [13] and HH re- 
spectively. These computed quantities are normalized by the shear velocity 
and compared with the DNS data [77]. It is noticed that the computed peak 
values of all the components of turbulent intensities are higher for Re* = 395 
than for Re* = 180. Moreover, the computed profiles are in good agreement 
with DNS. 

An important indication of how turbulence kinetic energy is transferred among 
the components is provided by the pressure-strain (PS) correlations. Their 
components are: PS X = (p d x u x ), PS y = (p d y u y ), and PS Z = (p d z u z ), where 
the prime denotes fluctuations and the brackets refer to averaging (along ho- 
mogeneous spatial directions and time). Figure [TBI shows the components of 
the computed PS correlations. They exhibit the expected behavior close to 
the wall, including the transfer of energy from the wall-normal component to 
the other two components near the wall - a phenomenon termed as splatting 
or impingement [TO] . 

It will be interesting to compare the computed results using the dynamic 
Smagorinsky SGS model with those using the "constant" Smagorinsky SGS 
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Fig. 13. Root-mean-square (rms) spanwise velocity fluctuations in fully-developed 
turbulent channel flow with a top slip surface at a shear Reynolds number of 
-Re* = 395 using a dynamic procedure with multiblock GLBE. Comparison with 
DNS of Moser, Kim and Mansour (1999). 
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Fig. 14. Root-mean-square(rms) wall normal velocity fluctuations in fully-devel- 
oped turbulent channel flow with a top slip surface at a shear Reynolds number of 
.Re* = 395 using a dynamic procedure with multiblock GLBE. Comparison with 
DNS of Moser, Kim and Mansour (1999). 
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Fig. 15. Components of pressure strain correlations in fully-developed turbulent 
channel flow with a top slip surface at a shear Reynolds number of Re* = 395 using 
a dynamic procedure with multiblock GLBE. 

model with an ad hoc van Driest damping function [9]. As an illustration, 
Fig. [16] shows a comparison of the components of turbulence intensities with 
both these SGS models. Also, plotted in symbols is the DNS data [77]. In the 
case of "constant" Smagorinsky SGS model, we consider C$ = 0.15 with the 



following damping function: A = 5 X 



exp 



where A + = 26. We find 



that even with the commonly used values of empirical constants, the results 
with the dynamical SGS model is in better agreement with the DNS. It may be 
noted that without the use of damping function, the "constant" Smagorinsky 
SGS model would have given grossly inaccurate results for this problem. In 
any case, the important point is that the dynamic SGS model, when used with 
the GLBE for LES, does not require any a priori information regarding the 
values or the behavior of the model coefficient, which are computed from the 
local behavior of the turbulent fields itself. 



6 Summary and Conclusions 



In this paper we discussed the incorporation of a dynamic procedure for large 
eddy simulation (LES) of complex turbulent flows using the generalized lat- 
tice Boltzmann equation (GLBE) with forcing term based on multiple relax- 
ation times. The use of dynamic procedure, which samples super-grid or "test" 
scale turbulent dynamics and assuming scale-invariance between two levels to 
compute the model coefficient, largely eliminates empiricism required to close 
the SGS turbulence models, particularly for inhomogeneous and anisotropic 
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Fig. 16. Comparison between constant (dashed) and dynamic (solid) Smagorin- 
sky models implemented in a multiblock GLBE framework for LES of fully-devel- 
oped turbulent channel flow with a top slip surface at a shear Reynolds number 
of -Re* = 395. Symbols are DNS data of Moser, Kim and Mansour (1999). For the 
constant Smagorinksy model, C s = 0.15 with van Driest damping using a coefficient 
of A + = 26 were considered. 

turbulent flows. Eddy-viscosity type dynamic Smagorinsky model (DSM) is 
employed as the SGS model in this work. The solution of the GLBE, which 
has significantly improved fidelity and numerical stability when compared to 
the SRT-LBE due to its ability to separate relaxation time scales of hydrody- 
namic and kinetic modes, yields the grid-filtered hydrodynamic fields, when 
"effective" relaxation times based on the dynamic SGS model is employed. 
Variable resolutions, to accommodate different turbulent length scales with 
the finest near the walls, are considered by the use of conservative, locally 
refined multiblock gridding approach. The grid-filtered strain-rate tensor used 
in the SGS model is computed locally from the non-equilibrium moments of 
the GLBE, which has superior stability characteristics when compared to the 
use of simple finite-differencing of the velocity field on coarser grids. To allow 
efficient implementation, discrete trapezoidal filters are employed to obtain 
the test-filtered quantities in the dynamic procedure. 

Numerical simulation results using the DSM are obtained for the canonical 
fully-developed turbulent channel flow at shear Reynolds numbers Re* of 180 
and 395. It is demonstrated that the use of the multiblock GLBE in con- 
junction with dynamic procedure for the SGS model is able to automatically 
compute the values and variation of the model coefficient C . In particular, 
without the use of any ad hoc approach, the computations reproduce the value 
of C equal to 0.024 (or Cs = 0.155) in the bulk flow region, which reduces 
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and approaches to zero towards to the wall, reflecting near-wall turbulence 
physics. The computed turbulent statistics such as the root-mean-square ve- 
locity fluctuations and Reynolds stresses for both cases of Reynolds numbers 
are in good agreement with prior direct numerical simulations (DNS) and ex- 
perimental data, as applicable. In particular, the peaks of the second-order 
statistics scale to higher values at higher .Re*, which the computations are 
found to reproduce quantitatively well. While the focus of this paper, as a 
first step, is on the use of dynamic procedure for eddy-viscosity based models 
in the LBM framework, it also outlines the details of implementation of more 
general scale-similarity based dynamic SGS models. It appears that the use of 
dynamic procedure in the LBM, particularly with the GLBE, is a promising 
approach for LES of complex turbulent flows. 
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A Moments, Equilibrium Moments, Moment-projections of Source 
Terms for D3Q19 Model 



The components of the various elements in the moments are as follows [56J: 

fo = P, fx = e, h = e 2 , fa = jx, h = Qx,fs = jy, /e = q v , h = jz, fs = q z ,fo = 

3p xx , flO — ^xx,fll — Pww, fl2 = fftuw, fl3 = Pxy, f 14 = Pyz, f 15 = Pxz, fvo = 

m x, fn = m y , fis — m z- Here, p is the density, e and e 2 represent kinetic 
energy that is independent of density and square of energy, respectively; j x , 
j y and j z are the components of the momentum, i.e. j x = pu x , j y = pu y , 
jz = pu z , q x , q y , q z are the components of the energy flux, and p xx , p xy , p yz 
and p xz are the components of the symmetric traceless viscous stress tensor; 
The other two normal components of the viscous stress tensor, p yy and p zz , 
can be constructed from p xx and p ww , where p ww = p yy — p zz . Other moments 
include ir xx , n ww , m x , m y and m z . The first two of these moments have the 
same symmetry as the diagonal part of the traceless viscous tensor pij, while 
the last three vectors are parts of a third rank tensor, with the symmetry of 

JkPmn- 

The corresponding components of the equilibrium moments, which are func- 
tions of conserved moments, i.e. density p and momentum fol- 
lows 1 561: 



fo q = p, f? = e eq = -np+i9^, f* = = 3 P -^^, if = j x , it = 

n eq — _2 • feq _ ■ ?eq eq _ 2 • 7eq _ ■ feq eq _ 2 • 7eq 

Hx — J5 — Jy, JQ — Hy — ^JyiJl — Jz, J8 — Hz ~ 3725/9 — 

o eq _ L x _ J _ J _ \ req _ o eq _ o(_l eq\ feq _ eq _ [h 3 Z \ f eq _ 
°Pxx — p iJlO — on xx — \ -iPxxJiJll — Pww — p )Jl2 — 

eq _ l n eq ?eq — n eq _ jxjy ?eq — n eq _ jyjz feq _ eq _ jxjz feq 
2P1VWIJI3 — Pxy — "Ti/14 — Pyz ~ ~~TiJ15 ~ Pxz ~ ~^~iJr 



71 ww — 2^w'wtJ13 — Vxy ~ p 1 J 14 — fyz ~ p ' J 15 — Pxz ~ p > J 16 

o,/r 7 9 = o,/s = o. 

The components of the source terms in moment space are functions of external 
force F and velocity fields ~Tt, respectively, as follows [62] : 



So = 0, S t = 38(F x u x + F R u y + F z u z ),S 2 = -ll(F x u x + F y u^ + F z u z ), S 3 = 
Fxi S^ = —^F x , Sn = F y , Sq = — ^F y , 5*7 = F z , Ss = —^F z , Sg = 2{2F x u x — 

FyUy ~ F Z U Z ),S W ^= ~(2F X U X - FyUy ~J^ Z U z ) , S U = 2(FyUy - F Z U z ) , 5 12 = 
-(FyUy^- F Z U Z \S 13 = (F X Uy + FyU X ),S 14 = (FyU Z + F Z Uy),S 15 = (F X U Z + 

F z u x ), S w = 0, #i7 = 0, Sis = 0. 

The components of the strain rate tensor used in subgrid scale (SGS) tur- 
bulence model at the grid-filter level can be written explicitly in terms of 
non-equilibrium moments augmented by moment-projections of source terms 
as 



28 



S r 



yy 



s z 



s 



xy ■ 



s„ 



yz 



S X z 

where 



1 



., Sl h { r 9) + 19^' 

Sop 1 
(bp L 
7bp L 



-(neq) 



19 (s 9 ht eq) - 3suhtr 9) ) 



(neq) 
11 



-— s 13 h 
2p 

3 



(neq) 
13 5 



2p 
3 

2p 



! '14 ) 



s 15 h 



(neq) 
15 5 



1 



t(neq) _ 7 _ Teq , rr 

a; —Jet Ja ' 2 



a G {1,9,11,13,14,15} 



(A.l) 
(A.2) 
(A.3) 
(A.4) 
(A.5) 
(A.6) 



(A.7) 



Here, / a , f* q and 5 Q are components of moments, their local equilibria, and 
moment-projections of source terms due to external forces, respectively, which 
are given above. s a are elements of the collision matrix A = diag(so, s\, . . . , sx$) 
in moment space. The expressions for the strain rate tensor are generalizations 
of those given in Yu et al. [50J 



B Incorporation of Dynamic Mixed Model (DMM) and Dynamic 
Two Parameter (DTM) Subgrid Scale Models in Lattice-Boltzmann 
Method 



It is well recognized that improvements to the dynamic Smagorinsky model 
(DSM) is required to properly account for various aspects of subgrid scale 
turbulence physics in general situations. For example, while DSM is generally 
adequate for performing large-eddy simulation (LES) of bulk as well as near- 
wall turbulent flows, it is less satisfactory near the free-slip or free- surfaces. 
This is presumably because of the assumption in the DSM that turbulence 
dissipation is in equilibrium with production, leading to the eddy-viscosity 
type models, which implies the alignment of the principal axes of the SGS 
stress tensor with those of the resolved strain rate tensor. Moreover, when 
the dynamic coefficient C is computed locally, it exhibits large fluctuations 
leading to excessive energy backscatter and numerical instability. Thus, the 
DSM can accurately predict mean flow quantities when an averaged value of 
the coefficient is employed, but is inadequate for representation of the local 
quantities. 



29 



Thus, Zang et al. (1993) [16] proposed to circumvent these problems by em- 
ploying a mixed model, which also contains contributions from the scale- 
similarity based stresses or the modified Leonard stress terms (as originally 
proposed by Bardina et al. (1980) [T5]) and do not assume the alignment of the 
SGS stress and the resolved strain rate tensors. The proportionality constant 
for the Leonard stress terms is set to 1.0 in this model. By alleviating the bur- 
den on the variation of the coefficient C to represent turbulence physics, this 
dynamic mixed model (DMM) provided improved results, with significantly 
less fluctuations in the model coefficient. A further improvement, particularly 
in the context of LES of free-surface turbulence, was provided by the dynamic 
two-parameter model (DTM) of Salvetti and Banerjee (1995) [17]. In this 
model, an additional parameter K is introduced as a proportionality constant 
to the modified Leonard stress term in the mixed based model. Both the co- 
efficients C and K can then be obtained locally by means of a least-squares 
technique, as formulated by Lilly (1992) [13] . We now discuss the develop- 
ment of a procedure for incorporation of DTM in the framework of lattice 
Boltzmann method (LBM). 

The anisotropic part of the SGS stress in the DTM can be represented as 
r ij - 5 -^T kk = -2Cf^\'3\'S ij + K 



j m 



6 ij j 



kk 



(B.l) 



where the first term on the right-hand-side (RHS) is due to the eddy-viscosity 
Smagorinsky model, while the second term is the modified Leonard stress term 
Ljj, which is given by 



(B.21 



The modified Leonard stress is a known quantity from the resolved field. At 
the test filter level, the SGS stress can be written as 
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Using Germano's identity (1991) [12], we obtain 
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which leads to 



3 



T %3 T 

J^in —J^kk 



-2CA Ma + K 



H a - -jH, 



kk 



(B.5) 



30 



where 



M ij = a?\130 ij -?tff ij . 



Here, all the above tensors, L y - , and M y -, are known from the resolved 

field and a = A/ A is the ratio of test-filter and grid-filter sizes. To obtain 
the constants C and K, the least squares technique in the context of dynamic 
SGS modeling as proposed by Lilly (1992) [TJ] is employed. Thus, the con- 
stants are obtained by minimizing the following function with respect to these 
parameters: 
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With dQ/dC = and dQ/dK = 0, we obtain 
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It may be noted that when DMM is employed, with K = 1, the Smagorinsky 
coefficient becomes 
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Once the local values of the coefficients C and K are known, the SGS stress can 
be incorporated within the framework of lattice Boltzmann method (LBM) as 
follows: The first part of the SGS stress in DTM provides the eddy viscosity 



v t = CA 2 \S\ 



(B.10) 



which can be added to the molecular viscosity u Q , obtained from the statement 
of the problem, to yield the total viscosity v , i.e., v — Vq + v t can then be 
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used to compute relaxation times (sio, S12, su, S15, sie) i n the MRT-LBE, as 
discussed in the main text of this paper for the DSM. The second part of the 
SGS stress can be introduced as forcing terms in the MRT-LBE. Thus, the 
modified Leonard stress can be written as the following Cartesian components 
of the effective force: 



■pmL 
TpmL 
jpmL 



d x (KL™)+d y (KL™)+d z (KLZ) 
d x (KL™)+d y (KL™)+d z (KL 
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The spatial derivatives needed in the forcing terms can be evaluated by em- 
ploying isotropic spatial discrete derivative operators based on lattice stencils. 

— imL — * 

Such forcing terms F , where they are added to F enter the MRT-LBE 
through the source terms in the moment space S a , defined at the end of Ap- 
pendix [Al 



It may be noted that the grid-filtered fields can be obtained directly from 
the solution of the MRT-LBE. On the other hand, the double-grid filtered 
<p as well as the test-filtered value 0, required in the DMM as well as the 
DSM, can be obtained from the repeated application of trapezoidal discrete 
filters for the grid and test volumes, respectively, successively for each spatial 
dimension [TO]. Thus, at grid node (i,j, k), for the double-grid filter: 



0(i,j,fc) = 1/8 + 60(ij,fc) + 4>(i-i,j,k)) > 

0(i,i,fc) = 1/ 8 (0(ij,fc+i) + 64>(i,j,k) + 0(ij,fc-i)) > 
and for the test-filter: 



(P(i,j,k) = V 4 (0(ij+i,k) + 2 ^( 4 ,i,fc) + <%j-i,fc) J > 

%,m = i/4 (% lfc +i) + 2$%j, k ) + • 

It may be noted that implementation of this approach in the MRT-LBE within 
the context of multiblock grids would require careful consideration of infor- 
mation exchange between different grid levels in the grid transition regions. 
Implementation and assessment of dynamic scale-similarity based SGS models 
in the GLBE framework are subjects of future investigations. 
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